!>
!>
!>     Purpose: This program evaluates the three-dimensional free-surface
!>              Green function and derivatives in deep water
!>
!>              Code Original Author: Hui Liang       created on  2016.12.26
!>
!>  License:
!>
!>
!>    This routine is a free software package: you can redistribute it and/or modify it
!>    under the terms of the GNU Lesser General Public License as published by the
!>    Free Software Foundation, either version 3 of the License, or (at your option)
!>    any later version.
!>
!>    You should have received a copy of the GNU General Public License (either V3
!>    or later), along with this routine. If not, see <http://www.gnu.org/licenses/>.
!>
!>  Modified on:
!>
!>    January 02, 2018
!>
!>  Reference:
!>
!>
!>        [1]        H. Wu, C. Zhang, Y. Zhu, W. Li, D. Wan, F. Noblesse,
!>                A global approximation to the Green function for
!>                diffraction radiation of water waves,
!>                Eur. J. Mech. B Fluids 65 (2017) 54-64.
!>
!>        [2]        H. Liang, H. Wu, F. Noblesse,
!>                Validation of a global approximation for
!>                wave diffraction-radiation in deep water,
!>                Appl. Ocean Res. 74 (2018) 80-86.
!>
!>   Remarks
!>
!>        The local-flow component is approximated by mean of the global
!>        approximations [1]. The computations reported in [2] provides
!>        strong evidence that the global approximations are sufficiently
!>        accurate to compute linear and second-order wave loads in practice.
!>
!>
!>        It should be noted that the Rankine source term -2/d appeared in L_z
!>        given by (8a) in Ref. [2] is not evaluated here (there is a typo in (8a)).
!>        These Rankine source components are required to evaluate in another routine
!>        due to strong singular behavior.
!>
!>        For any questions, please contact: lianghuistar@gmail.com
!>
!>   We define the flow-field point (x,y,z) and source point (xi,eta,zeta).
!>   Please note that all variables are non-dimensionalized with respect to
!>        the wavenumber k0
!>
!>   The Green function is defiend as
!>
!>   \( G = -1/r-1/d+GF \),
!>
!>   where
!>
!>   \( r = sqrt((x-xi)^2+(y-eta)^2+(z-zeta)^2) \); <br>
!>   \( d = sqrt((x-xi)^2+(y-eta)^2+(z+zeta)^2) \).
!>
!>   Parameters:
!>
!>       Input:      dx  --- x-xi
!>                   dy  --- y-eta
!>                   vv  --- z+zeta
!>
!>       Output:     GF(0)  --- free surface itself \( GF \)
!>                   GF(1)  --- x-derivative \( GF_x \)
!>                   GF(2)  --- y-derivative \( GF_y \)
!>                   GF(3)  --- part of z-derivative \( GF_z = GF(3) - 2/d \)
module infin_green_3d_m

    use, intrinsic :: iso_fortran_env, only: rk => real64
    implicit none
    private

    public :: HavelockGF

    real(rk) :: pi = acos(-1.0_rk)
    real(rk) :: gama = 0.5772156649_rk
    complex(rk) :: Im = (0.0_rk, 1.0_rk)

    real(rk) :: GF_alpha
    real(rk) :: GF_beta
    real(rk) :: GF_sigma
    real(rk) :: GF_rho
    real(rk) :: GF_dd

contains

    !> This program evaluates the three-dimensional free-surface
    !> Green function and derivatives in deep water <br>
    !> 该程序评估深水中的三维自由表面格林函数及其导数
    impure subroutine HavelockGF(dx, dy, vv, GF)
        real(rk), intent(in) :: dx          !! \( x-xi \)
        real(rk), intent(in) :: dy          !! \( y-eta \)
        real(rk), intent(in) :: vv          !! \( z+zeta \)
        complex(rk), intent(out) :: GF(0:3) !! free surface itself \( GF \) <br>
                                            !! GF(0)  --- free surface itself \( GF \) <br>
                                            !! GF(1)  --- x-derivative \( GF_x \) <br>
                                            !! GF(2)  --- y-derivative \( GF_y \) <br>
                                            !! GF(3)  --- part of z-derivative \( GF_z = GF(3) - 2/d \)

        real(rk) :: hh
        complex(rk) :: GFh

        hh = dsqrt(dx*dx + dy*dy)
        call GF_DivParameters(hh, vv)

        GF(0) = GF_Func_L0(hh, vv) + GF_Func_W(hh, vv)
        GFh = GF_Func_Ls(hh, vv) + GF_Func_Wh(hh, vv)

        if (hh > 1.0d-6) then
            GF(1) = GFh*dx/hh
            GF(2) = GFh*dy/hh
        else
            GF(1) = 0.0_rk
            GF(2) = 0.0_rk
        end if

        GF(3) = GF(0)

    end subroutine HavelockGF

    impure subroutine GF_DivParameters(hh, vv)
        real(rk), intent(in) :: hh, vv

        GF_dd = dsqrt(hh*hh + vv*vv)
        GF_alpha = -vv/GF_dd
        GF_beta = hh/GF_dd
        GF_sigma = hh/(GF_dd - vv)
        GF_rho = GF_dd/(1.0_rk + GF_dd)

    end subroutine GF_DivParameters

    pure function GF_Func_L0(hh, vv)
        real(rk), intent(in) :: hh, vv
        integer :: ii
        real(rk) :: PP, Lp, Cal_L0
        real(rk) :: GF_Func_L0

        PP = dlog(0.5_rk*(GF_dd - vv)) + gama - 2.0_rk*GF_dd*GF_dd
        PP = dexp(vv)*PP
        PP = PP + GF_dd*GF_dd - vv

        Lp = GF_Func_Lp(hh, vv)

        GF_Func_L0 = 2.0_rk*PP/(1.0_rk + GF_dd**3) + 2.0_rk*Lp

    end function GF_Func_L0

    pure function GF_Func_Lp(hh, vv)
        real(rk), intent(in) :: hh, vv
        real(rk) :: A, B, C, D, RR, Lp
        real(rk) :: GF_Func_Lp

        A = GF_FuncA(GF_rho)
        B = GF_FuncB(GF_rho)
        C = GF_FuncC(GF_rho)
        D = GF_FuncD(GF_rho)

        RR = (1.0_rk - GF_beta)*A
        RR = RR - GF_beta*B
        RR = RR - GF_alpha*C/(1.0_rk + 6.0_rk*GF_alpha*GF_rho*(1.0_rk - GF_rho))
        RR = RR + GF_beta*(1.0_rk - GF_beta)*D

        GF_Func_Lp = GF_rho*(1.0_rk - GF_rho)**3*RR

    end function GF_Func_Lp

    pure function GF_Func_W(hh, vv)
        real(rk), intent(in) :: hh, vv
        real(rk) :: H0, J0
        complex(rk) :: GF_Func_W

        H0 = StruveH0(hh)
        J0 = BesselJ0(hh)

        GF_Func_W = 2.0_rk*pi*(H0 - Im*J0)*dexp(vv)

    end function GF_Func_W

    pure function GF_FuncA(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_FuncA
        
        real(rk), parameter :: A0 = 1.21_rk
        real(rk), parameter :: A1 = -13.328_rk
        real(rk), parameter :: A2 = +215.896_rk
        real(rk), parameter :: A3 = -1763.96_rk
        real(rk), parameter :: A4 = +8418.94_rk
        real(rk), parameter :: A5 = -24314.21_rk
        real(rk), parameter :: A6 = +42002.57_rk
        real(rk), parameter :: A7 = -41592.9_rk
        real(rk), parameter :: A8 = 21859.0_rk
        real(rk), parameter :: A9 = -4838.6_rk

        GF_FuncA = (A8 + A9*tt)*tt
        GF_FuncA = (A7 + GF_FuncA)*tt
        GF_FuncA = (A6 + GF_FuncA)*tt
        GF_FuncA = (A5 + GF_FuncA)*tt
        GF_FuncA = (A4 + GF_FuncA)*tt
        GF_FuncA = (A3 + GF_FuncA)*tt
        GF_FuncA = (A2 + GF_FuncA)*tt
        GF_FuncA = (A1 + GF_FuncA)*tt
        GF_FuncA = A0 + GF_FuncA

    end function GF_FuncA

    pure function GF_FuncB(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_FuncB
        
        real(rk), parameter :: B0 = 0.938_rk
        real(rk), parameter :: B1 = +5.737_rk
        real(rk), parameter :: B2 = -67.92_rk
        real(rk), parameter :: B3 = +796.534_rk
        real(rk), parameter :: B4 = -4780.77_rk
        real(rk), parameter :: B5 = +17137.74_rk
        real(rk), parameter :: B6 = -36618.81_rk
        real(rk), parameter :: B7 = +44894.06_rk
        real(rk), parameter :: B8 = -29030.24_rk
        real(rk), parameter :: B9 = +7671.22_rk

        GF_FuncB = (B8 + B9*tt)*tt
        GF_FuncB = (B7 + GF_FuncB)*tt
        GF_FuncB = (B6 + GF_FuncB)*tt
        GF_FuncB = (B5 + GF_FuncB)*tt
        GF_FuncB = (B4 + GF_FuncB)*tt
        GF_FuncB = (B3 + GF_FuncB)*tt
        GF_FuncB = (B2 + GF_FuncB)*tt
        GF_FuncB = (B1 + GF_FuncB)*tt
        GF_FuncB = B0 + GF_FuncB

    end function GF_FuncB

    pure function GF_FuncC(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_FuncC

        real(rk), parameter :: C0 = 1.268_rk
        real(rk), parameter :: C1 = -9.747_rk
        real(rk), parameter :: C2 = +209.653_rk
        real(rk), parameter :: C3 = -1397.89_rk
        real(rk), parameter :: C4 = +5155.67_rk
        real(rk), parameter :: C5 = -9844.35_rk
        real(rk), parameter :: C6 = +9136.4_rk
        real(rk), parameter :: C7 = -3272.62_rk
        

        GF_FuncC = (C6 + C7*tt)*tt
        GF_FuncC = (C5 + GF_FuncC)*tt
        GF_FuncC = (C4 + GF_FuncC)*tt
        GF_FuncC = (C3 + GF_FuncC)*tt
        GF_FuncC = (C2 + GF_FuncC)*tt
        GF_FuncC = (C1 + GF_FuncC)*tt
        GF_FuncC = C0 + GF_FuncC

    end function GF_FuncC

    pure function GF_FuncD(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_FuncD
            
        real(rk), parameter :: D0 = 0.632_rk
        real(rk), parameter :: D1 = -40.97_rk
        real(rk), parameter :: D2 = 667.16_rk
        real(rk), parameter :: D3 = -6072.07_rk
        real(rk), parameter :: D4 = 31127.39_rk
        real(rk), parameter :: D5 = -96293.05_rk
        real(rk), parameter :: D6 = 181856.75_rk
        real(rk), parameter :: D7 = -205690.43_rk
        real(rk), parameter :: D8 = 128170.2_rk
        real(rk), parameter :: D9 = -33744.6_rk
        

        GF_FuncD = (D8 + D9*tt)*tt
        GF_FuncD = (D7 + GF_FuncD)*tt
        GF_FuncD = (D6 + GF_FuncD)*tt
        GF_FuncD = (D5 + GF_FuncD)*tt
        GF_FuncD = (D4 + GF_FuncD)*tt
        GF_FuncD = (D3 + GF_FuncD)*tt
        GF_FuncD = (D2 + GF_FuncD)*tt
        GF_FuncD = (D1 + GF_FuncD)*tt
        GF_FuncD = D0 + GF_FuncD

    end function GF_FuncD

    pure function GF_Func_Ls(hh, vv)
        real(rk), intent(in) :: hh, vv
        real(rk) :: PS, QS, Lsp
        real(rk) :: GF_Func_Ls

        PS = (GF_beta + hh)/(GF_dd - vv)
        PS = PS - 2.0_rk*GF_beta + 2.0_rk*dexp(vv)*GF_dd - hh

        QS = dexp(-GF_dd)*(1.0_rk - GF_beta)
        QS = QS*(1.0_rk + GF_dd/(1.0_rk + GF_dd**3))

        Lsp = GF_Func_Lsp(hh, vv)

        GF_Func_Ls = 2.0_rk*PS/(1.0_rk + GF_dd**3) - 4.0_rk*QS + 2.0_rk*Lsp

    end function GF_Func_Ls

    pure function GF_Func_Lsp(hh, vv)
        real(rk), intent(in) :: hh, vv
        real(rk) :: A, B, C, RR
        real(rk) :: GF_Func_Lsp

        A = GF_dFuncA(GF_rho)
        B = GF_dFuncB(GF_rho)
        C = GF_dFuncC(GF_rho)

        RR = GF_beta*A
        RR = RR - (1.0_rk - GF_alpha)*B
        RR = RR + GF_beta*(1.0_rk - GF_beta)*GF_rho*(1.0_rk - 2.0_rk*GF_rho)*C

        GF_Func_Lsp = GF_rho*(1.0_rk - GF_rho)**3*RR

    end function GF_Func_Lsp

    pure function GF_Func_Wh(hh, vv)
        real(rk), intent(in) :: hh, vv
        real(rk) :: H1, J1
        complex(rk) :: GF_Func_Wh

        H1 = StruveH1(hh)
        J1 = BesselJ1(hh)

        GF_Func_Wh = 2.0_rk*pi*(2.0_rk/pi - H1 + Im*J1)*dexp(vv)

    end function GF_Func_Wh

    pure function GF_dFuncA(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_dFuncA
        
        real(rk), parameter :: A0 = 2.948_rk
        real(rk), parameter :: A1 = -24.53_rk
        real(rk), parameter :: A2 = +249.69_rk
        real(rk), parameter :: A3 = -754.85_rk
        real(rk), parameter :: A4 = -1187.71_rk
        real(rk), parameter :: A5 = +16370.75_rk
        real(rk), parameter :: A6 = -48811.41_rk
        real(rk), parameter :: A7 = +68220.87_rk
        real(rk), parameter :: A8 = -46688.0_rk
        real(rk), parameter :: A9 = +12622.25_rk
        

        GF_dFuncA = (A8 + A9*tt)*tt
        GF_dFuncA = (A7 + GF_dFuncA)*tt
        GF_dFuncA = (A6 + GF_dFuncA)*tt
        GF_dFuncA = (A5 + GF_dFuncA)*tt
        GF_dFuncA = (A4 + GF_dFuncA)*tt
        GF_dFuncA = (A3 + GF_dFuncA)*tt
        GF_dFuncA = (A2 + GF_dFuncA)*tt
        GF_dFuncA = (A1 + GF_dFuncA)*tt
        GF_dFuncA = A0 + GF_dFuncA

    end function GF_dFuncA

    pure function GF_dFuncB(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_dFuncB
        
        real(rk), parameter :: B0 = 1.11_rk
        real(rk), parameter :: B1 = 2.894_rk
        real(rk), parameter :: B2 = -76.765_rk
        real(rk), parameter :: B3 = 1565.35_rk
        real(rk), parameter :: B4 = -11336.19_rk
        real(rk), parameter :: B5 = 44270.15_rk
        real(rk), parameter :: B6 = -97014.11_rk
        real(rk), parameter :: B7 = 118879.26_rk
        real(rk), parameter :: B8 = -76209.82_rk
        real(rk), parameter :: B9 = 19923.28_rk

        GF_dFuncB = (B8 + B9*tt)*tt
        GF_dFuncB = (B7 + GF_dFuncB)*tt
        GF_dFuncB = (B6 + GF_dFuncB)*tt
        GF_dFuncB = (B5 + GF_dFuncB)*tt
        GF_dFuncB = (B4 + GF_dFuncB)*tt
        GF_dFuncB = (B3 + GF_dFuncB)*tt
        GF_dFuncB = (B2 + GF_dFuncB)*tt
        GF_dFuncB = (B1 + GF_dFuncB)*tt
        GF_dFuncB = B0 + GF_dFuncB

    end function GF_dFuncB

    pure function GF_dFuncC(tt)
        real(rk), intent(in) :: tt
        real(rk) :: GF_dFuncC
        
        real(rk), parameter :: C0 = 14.19_rk
        real(rk), parameter :: C1 = -148.24_rk
        real(rk), parameter :: C2 = 847.8_rk
        real(rk), parameter :: C3 = -2318.58_rk
        real(rk), parameter :: C4 = 3168.35_rk
        real(rk), parameter :: C5 = -1590.27_rk
        

        GF_dFuncC = (C4 + C5*tt)*tt
        GF_dFuncC = (C3 + GF_dFuncC)*tt
        GF_dFuncC = (C2 + GF_dFuncC)*tt
        GF_dFuncC = (C1 + GF_dFuncC)*tt
        GF_dFuncC = C0 + GF_dFuncC

    end function GF_dFuncC

    pure function BesselJ0(xx)
        real(rk), intent(in) :: xx
        real(rk) :: yy, y2, f0, theta0
        real(rk) :: BesselJ0
            
        real(rk), parameter :: P0 = 0.999999999_rk
        real(rk), parameter :: P1 = -2.249999879_rk
        real(rk), parameter :: P2 = 1.265623060_rk
        real(rk), parameter :: P3 = -0.316394552_rk
        real(rk), parameter :: P4 = 0.044460948_rk
        real(rk), parameter :: P5 = -0.003954479_rk
        real(rk), parameter :: P6 = 0.000212950_rk
        real(rk), parameter :: R0 = 0.79788454_rk
        real(rk), parameter :: R1 = -0.00553897_rk
        real(rk), parameter :: R2 = 0.00099336_rk
        real(rk), parameter :: R3 = -0.00044346_rk
        real(rk), parameter :: R4 = 0.00020445_rk
        real(rk), parameter :: R5 = -0.00004959_rk
        real(rk), parameter :: S1 = -0.04166592_rk
        real(rk), parameter :: S2 = 0.00239399_rk
        real(rk), parameter :: S3 = -0.00073984_rk
        real(rk), parameter :: S4 = 0.00031099_rk
        real(rk), parameter :: S5 = -0.00007605_rk
        

        if (xx <= 3.0_rk) then
            yy = (xx/3.0_rk)**2
            BesselJ0 = P0 + (P1 + (P2 + (P3 + (P4 + (P5 + P6*yy)*yy)*yy)*yy)*yy)*yy
        else
            yy = 3.0_rk/xx
            y2 = yy**2

            f0 = R0 + (R1 + (R2 + (R3 + (R4 + R5*y2)*y2)*y2)*y2)*y2

            theta0 = xx - 0.25_rk*pi &
                     + (S1 + (S2 + (S3 + (S4 + S5*y2)*y2)*y2)*y2)*yy

            BesselJ0 = f0*dcos(theta0)/dsqrt(xx)
        end if

    end function BesselJ0

    pure function BesselJ1(xx)
        real(rk), intent(in) :: xx
        real(rk) :: yy, y2, f1, theta1
        real(rk) :: BesselJ1
            
        real(rk), parameter :: P0 = 0.500000000_rk
        real(rk), parameter :: P1 = -0.562499992_rk
        real(rk), parameter :: P2 = 0.210937377_rk
        real(rk), parameter :: P3 = -0.039550040_rk
        real(rk), parameter :: P4 = 0.004447331_rk
        real(rk), parameter :: P5 = -0.000330547_rk
        real(rk), parameter :: P6 = 0.000015525_rk
        real(rk), parameter :: R0 = 0.79788459_rk
        real(rk), parameter :: R1 = 0.01662008_rk
        real(rk), parameter :: R2 = -0.00187002_rk
        real(rk), parameter :: R3 = 0.00068519_rk
        real(rk), parameter :: R4 = -0.00029440_rk
        real(rk), parameter :: R5 = 0.00006952_rk
        real(rk), parameter :: S1 = 0.12499895_rk
        real(rk), parameter :: S2 = -0.00605240_rk
        real(rk), parameter :: S3 = 0.00135825_rk
        real(rk), parameter :: S4 = -0.00049616_rk
        real(rk), parameter :: S5 = 0.00011531_rk
        

        if (xx <= 3.0_rk) then
            yy = xx/3.0_rk
            y2 = yy*yy
            BesselJ1 = P0 + (P1 + (P2 + (P3 + (P4 + (P5 + P6*y2)*y2)*y2)*y2)*y2)*y2

            BesselJ1 = BesselJ1*xx
        else
            yy = 3.0_rk/xx
            y2 = yy*yy
            f1 = R0 + (R1 + (R2 + (R3 + (R4 + R5*y2)*y2)*y2)*y2)*y2

            theta1 = xx - 0.75_rk*pi &
                     + (S1 + (S2 + (S3 + (S4 + S5*y2)*y2)*y2)*y2)*yy

            BesselJ1 = f1*dcos(theta1)/dsqrt(xx)
        end if

    end function BesselJ1

    pure function BesselY0(xx)
        real(rk), intent(in) :: xx
        real(rk) :: f0, theta0
        real(rk) :: BesselY0

        if (xx <= 3.0_rk) then
            BesselY0 = (2.0_rk/pi)*dlog(xx/2.0_rk)*BesselJ0(xx) &
                       + 0.367466907_rk &
                       + 0.605593797_rk*(xx/3.0_rk)**2 &
                       - 0.743505078_rk*(xx/3.0_rk)**4 &
                       + 0.253005481_rk*(xx/3.0_rk)**6 &
                       - 0.042619616_rk*(xx/3.0_rk)**8 &
                       + 0.004285691_rk*(xx/3.0_rk)**10 &
                       - 0.000250716_rk*(xx/3.0_rk)**12
        else
            f0 = 0.79788454_rk &
                 - 0.00553897_rk*(3.0_rk/xx)**2 &
                 + 0.00099336_rk*(3.0_rk/xx)**4 &
                 - 0.00044346_rk*(3.0_rk/xx)**6 &
                 + 0.00020445_rk*(3.0_rk/xx)**8 &
                 - 0.00004959_rk*(3.0_rk/xx)**10

            theta0 = xx - 0.25_rk*pi &
                     - 0.04166592_rk*(3.0_rk/xx)**1 &
                     + 0.00239399_rk*(3.0_rk/xx)**3 &
                     - 0.00073984_rk*(3.0_rk/xx)**5 &
                     + 0.00031099_rk*(3.0_rk/xx)**7 &
                     - 0.00007605_rk*(3.0_rk/xx)**9

            BesselY0 = f0*dsin(theta0)/dsqrt(xx)

        end if

    end function BesselY0

    pure function BesselY1(xx)
        real(rk), intent(in) :: xx
        real(rk) :: f1, theta1
        real(rk) :: BesselY1

        if (xx <= 3.0_rk) then
            BesselY1 = (2.0_rk/pi)*(dlog(xx/2.0_rk)*BesselJ1(xx) - 1.0_rk/xx) &
                       + 0.07373571_rk*(xx/3.0_rk)**1 &
                       + 0.72276433_rk*(xx/3.0_rk)**3 &
                       - 0.43885620_rk*(xx/3.0_rk)**5 &
                       + 0.10418264_rk*(xx/3.0_rk)**7 &
                       - 0.01340825_rk*(xx/3.0_rk)**9 &
                       + 0.00094249_rk*(xx/3.0_rk)**11
        else
            f1 = 0.79788459_rk &
                 + 0.01662008_rk*(3.0_rk/xx)**2 &
                 - 0.00187002_rk*(3.0_rk/xx)**4 &
                 + 0.00068519_rk*(3.0_rk/xx)**6 &
                 - 0.00029440_rk*(3.0_rk/xx)**8 &
                 + 0.00006952_rk*(3.0_rk/xx)**10

            theta1 = xx - 3.0_rk*pi/4.0_rk &
                     + 0.12499895_rk*(3.0_rk/xx)**1 &
                     - 0.00605240_rk*(3.0_rk/xx)**3 &
                     + 0.00135825_rk*(3.0_rk/xx)**5 &
                     - 0.00049616_rk*(3.0_rk/xx)**7 &
                     + 0.00011531_rk*(3.0_rk/xx)**9

            BesselY1 = f1*dsin(theta1)/dsqrt(xx)

        end if

    end function BesselY1

    pure function StruveH0(xx)
        real(rk), intent(in) :: xx

        integer :: ii
        real(rk) :: P0, P1, P2
        real(rk) :: P3, P4, P5
        real(rk) :: a0, a1, a2, a3
        real(rk) :: b1, b2, b3
        real(rk) :: c1, c2
        real(rk) :: yy, StruveH0

        if (xx <= 3.0_rk) then

            yy = (xx/3.0_rk)**2

            P0 = +1.909859164_rk
            P1 = -1.909855001_rk
            P2 = +0.687514637_rk
            P3 = -0.126164557_rk
            P4 = +0.013828813_rk
            P5 = -0.000876918_rk

            StruveH0 = P0 + (P1 + (P2 + (P3 + (P4 + P5*yy)*yy)*yy)*yy)*yy

            StruveH0 = StruveH0*(xx/3.0_rk)

        else

            yy = (3.0_rk/xx)**2

            a0 = 0.99999906_rk
            a1 = 4.77228920_rk
            a2 = 3.85542044_rk
            a3 = 0.32303607_rk

            b1 = 4.88331068_rk
            b2 = 4.28957333_rk
            b3 = 0.52120508_rk

            c1 = 2.0_rk*(a0 + (a1 + (a2 + a3*yy)*yy)*yy)
            c2 = pi*xx*(1.0_rk + (b1 + (b2 + b3*yy)*yy)*yy)

            StruveH0 = c1/c2 + BesselY0(xx)

        end if

    end function StruveH0

    pure function StruveH1(xx)
        real(rk), intent(in) :: xx

        integer :: ii
        real(rk) :: P1, P2, P3
        real(rk) :: P4, P5, P6
        real(rk) :: a0, a1, a2, a3
        real(rk) :: b1, b2, b3
        real(rk) :: c1, c2, yy
        real(rk) :: StruveH1

        if (xx <= 3.0_rk) then

            yy = (xx/3.0_rk)**2

            P1 = +1.909859286_rk
            P2 = -1.145914713_rk
            P3 = +0.294656958_rk
            P4 = -0.042070508_rk
            P5 = +0.003785727_rk
            P6 = -0.000207183_rk

            StruveH1 = (P1 + (P2 + (P3 + (P4 + (P5 + P6*yy)*yy)*yy)*yy)*yy)*yy

        else

            yy = (3.0_rk/xx)**2

            a0 = 1.00000004_rk
            a1 = 3.92205313_rk
            a2 = 2.64893033_rk
            a3 = 0.27450895_rk

            b1 = 3.81095112_rk
            b2 = 2.26216956_rk
            b3 = 0.10885141_rk

            c1 = 2.0_rk*(a0 + (a1 + (a2 + a3*yy)*yy)*yy)
            c2 = pi*(1.0_rk + (b1 + (b2 + b3*yy)*yy)*yy)

            StruveH1 = c1/c2 + BesselY1(xx)

        end if

    end function StruveH1

end module infin_green_3d_m
